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ABSTRACT 

We perform a series of two-dimensional magnetohydrodynamic core-collapse 
simulations of rapidly rotating and strongly magnetized massive stars. To study 
the properties of magnetic explosions for a longer time stretch of postbounce 
evolution, we develop a new code under the framework of special relativity in- 
cluding a realistic equation of state with a multiflavor neutrino leakage scheme. 
Our results show the generation of the magnetically-dominated jets in the two 
ways. One is launched just after the core-bounce in a prompt way and another 
is launched at ~ 100 ms after the stall of the prompt shock. We find that the 
shock-revival occurs when the magnetic pressure becomes strong, due to the field 
wrapping, enough to overwhelm the ram pressure of the accreting matter. The 
critical toroidal magnetic fields for the magnetic shock-revival are found to be 
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universal of ~ 10 15 G behind the jets. We point out that the time difference be- 
fore the shock-revival has a strong correlation with the explosions energies. Our 
results suggest that the magnetically dominated jets are accompanied by the for- 
mation of the magnetars. Since the jets are mildly relativistic, we speculate that 
they might be the origin of some observed X-ray flashes. 

Subject headings: supernovae: collapse, rotation — magnetars: pulsars, magnetic 
field — methods: numerical — MHD — special relativity — gamma rays: bursts 



1. Introduction 



There have been growing evidences shedding lights on the relations between the high 
energy astrophysical phenomena and their origins. A number of host galaxies of long- 
duration gamma-ray bursts (GRBs) are recently identified as metal poor galaxies whose 



meta l icities are lower th an that of average massive star-forming galaxies (jSavaglio et al. 



20061 : IStanek et al.ll2006l . and reference therein). The preponderance of short-lived massive 
star formation in such young galaxies, as well as the identification of SN Ib/c light curves 
peaking after the bursts in a fe w cases, has provided strong suppor t for a massive stell ar col- 
lapse origin of the long GRBs (|Paczynskilll998l ; iGalama et al.lll998l ; IStanek et al.ll2003l ). The 
duration of the l ong GRBs may c orrespond to the accretion of debris falling into the central 
black hole (BH) (jPiro et al.lll998l ). which suggests the observational consequence of the BH 
formation likewise the supernova of neutron star formation. There is also a growing observa- 
tional evidence of super magnetized neutron stars with the magnetic fields of ~ 10 14 — 10 15 G, 



the so-called magnetars (IDuncan fc Thompson! (119921 ). see lLattimer fc Prakashl (120071 ) for a 
recent review) . The magnetic fields are determined by the measured period and derivative of 
period un der the assumption that t he spin-down is cause d due to the usual magnetic dipole 
radiation (jZhang k, Harding! I20QCH : lHarding fc Laill2006l ). Tentative detections of spectral 
features duri ng the burst phase also indicate B ~_ 10 15 G when inte rpreted as proton cy- 



clotron lines (jGavriil et al 



2002 : llbrahim et al.ll2003l : iRea et al.ll2003f ). Recently X-ray flash 



(XRF), which is a low energy analogue of t he GRB, is receiving great attentions as a pos- 
sible relevance to the magnetar formations (jMazzali et al.ll2006l : iToma et al.l 120071 ). A large 
amount of neutron rich Ni ejected by SN2006aj associated with XRF06 0218 is interpreted t o 
be the formation of such objects, not the black hole after the explosion (IMaeda et al.ll2007al ). 



So far a number of numerical simulations have been done towards the understanding of 
the formation mechanisms of these compact objects such as neutron stars, magnetars, and 
the black holes in combination with their possible consequences like GRBs and XRFs. The 
leading model for the long-duration GRBs is the collapsar model ( IMacFadyen fc Woosley 
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19991 ). In the model, the core of massive stars with significant angular momentum collapses 
into a black hole (BH). The neutrinos emitted from the rotation-supported accretion disk 
around the BH heat the matter of the funnel region of the disk, to launch the GRB outflows. 
The relativistic flows are expected t o ultimately form a fireball, which is good for the explana- 
tion of the observed afterglow (e.g.. |Piranl (1l999f ) ) . In addition, it is suggested that the strong 
magnetic fields in the cores of order of 10 15 G play also an active role both for driving the 
magn eto-dr i ven je t s and for extrac t ing a si gnificant amount of ener g y from the central engine 
(e.g., Illsoj Jl992h : lwheeler et all (boooh : iThompson et all tooj ): lUzdenskv fc MacFadven 
(j2007al ) and see references therein). 



In order to understand such scenarios, the ultimate necessity of the stellar core-collapse 
simulations is to perform the simulations tracing all the phases in a consistent manner 
starting from the stellar core-collapse, core-bounce, shock-stall, stellar explosion (phase 1) 
or BH formation and the formation of accretion disk (phase 2), energy deposition to the 
funnel region by neutrinos and/or magnetic fields (phase 3), to the launching of the fireballs 
(phase 4). Here for convenience we call each stage as phase 1, 2, etc. The requirement 
for the numerical modeling to this end is highly computationally expensive, namely the 
multidimensional MHD simulations not only with general relativity for handling the BH 
formation, but also with the multi-angle neutrino transfer for treating highly anisotropic 
neutrino radiation from the disks. So various approximative approaches to each phase have 
been undertaken. As we mention below, these studies are complimentary in the sense that 
the different epochs are focused on, with the different initial conditions for the numerical 
modeling being taken. 

In ad dition to the e l abora t e studies in the con ventional supernova context (see recent 



reviews for lKotake et al.l (120061 ) ; I Janka et al.l (120071 )). much attention has been paid recently 



and its possi 


o\e app 


ication to t 


re collapsars (Yamada & Sawai 12004; Takiwaki et al. 2004 




Kotake et al. 


2004bl: 


Sawai e 


t al. 


2005:lMat 


; et al. 20061; Moiseenko et al.l/ 


!006l; Obergaulinger et al. 


2006; 


Nishimura et al.ll2006; 


3urrows et al. 


2007; 


Cerda-Duran et al.l 2007 


; Scheidegger et al 




2007: 


Komissarov & Barkov 


2007). After the failed or weak explosion, the accretion to the 



central objects may lead to the formation of a BH (phase 2). Several general relativistic 
studies are on the line for the understanding of the hydrodynamics at the epoch of the BH 
formation, in which more massi ve progenitors (>~ 2 5M^) than those of the st udy in the 



phase 1 are generally employed (jShibata et al. 



20061 : ISekiguchi & Shibatall2007f ). Treating 



the BH as an absorbing boundary or using the fixed metric approaches, the numerical stud- 
ies of the phase 3 are concerned with the initiation of the outflows from the funnel region 
of the disk to the acceleration of the jets as a re sult of the neutrino heating and/or MHD 
processes till the jets become mildly relativistic (IKoide et al.l 11998c iMacFadyen fc Woosley 
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1999; 


Proga et a 


2003 


Nishikawa et al. 


2C 


05; De Villiers et al. 


2005: Krolik et al. 


2005 




Hawl< 


3V & Krolik 


2006; 


Vlizuno e 


t al. 200 




Fuiimoto et a . 


2006 


Uzdenskv & 


VlacFadver 


2006: 


McKinnev & Naravanl2007: 


Komissarov & McKinnev 


2007; 


Magataki et al. 


2007 


Suwa et al. 


20071 


) al; Barkov & Komissarovl2008). Numerical studies of the phase 4 are mainly concerned 



with the dynamics later on, namely, the jet propagation to the break out from the star, when 



the acceleration of the jets to the h igh Lorentz factor is expected (IStone fc Hardee 



Alov et al. 


2000; 


Z 


lang et al. 


2006; 


Mizuno et al. 


2007 





lang et al.l 120031 ; iLeismann et al.l 120051 ; iMizuta et al. 



200 



6; 



2000 



McKinnev 



Our previous study was devoted to the phase 1, in which we performed a series of 2D 
core-collapse simulations of rotating an d magnetized massive stars under the framework of 
the Newtonian magnetohydrodynamics (ITakiwaki et al.ll2004l ). We found that the magneto- 
driven jet-like shocks were launched from the protoneutron stars just after core-bounce. 
However at the moment, we were unable to follow the dynamics much later on until when the 
collimated jets reach further out from the center. The Alfven velocity of the jet propagating 
into the outer layer of the iron core can be estimated by the following simple order-of- 
magnitude estimation, 



B 



Va 



10 iU cm/s 



£/10 13 G 



cur 



with p and B being the typical density and magnetic field there. It can be readily inferred 
that the Alfven velocity can exceed the speed of light unphysically in the Newtonian sim- 
ulation. To avoid this problem we construct a new code under the framework of special 
relativistic magnetohydrodynamics. We take a wider parametric range for the strength of 
the rotation than that of our previous work. By so doing, we hope to study more systemati- 
cally than before how the strong magnetic fields and the rapid rotation affect the properties 
of the magnetic explosions. 

We summarize the numerical methods in section [2J Section 3 is devoted to the initial 
models. In section 4, we show the numerical results. In section 5, we summarize our study 
and discuss the implications of our model for the magnetars and the X-ray flashes. Details 
of the numerical scheme and the code tests are given in the appendix. 



2. Numerical Methods 

The results presented in this paper are calculated by the newly developed special rel- 
ativistic magnetohydrodynamic (SRMHD) code. The novel point of this code is that the 
detailed microphysical processes relevant for the stellar-core-collapse simulations are also 
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coupled to the magneto-hydrodynamics (MHD). We briefly summarize the numerical meth- 
ods in the following. 



The MHD part of the code is based on the formalism of iDe Villiers et al.l (120031 ). Before 
going to the basic equations, we write down the definition of the primary code variables. 
The state of the relativistic fluid element at each point in the space time is described by its 
density, p; specific energy, e; velocity, v % \ and pressure, p. And the magnetic field in the rest 
frame of the fluid is described by the 4-vector \f^nb^ = *F tiy U y) where *F^ U is the dual of 



the electro-magnetic field strength tensor and U u is the 4-velocity. 

After some mathematical procedures presented in Appendix |A] the basic equations of 
SRMHD are described as follows: 

f + j=e lV jDv' = (2) 

dE 1 dW v 

^ + -^vW = -P^-^W-L, (3) 

dSi m b% +^d j ^(sy-b i v) = -\{ P h{Wv k f-{b k f)d a kk 

- [phW 2 - 6* 2 ) di<S> 

-9i(p + ^P\ (4) 



J +dj{Wv j b l -WvV) = (5) 



dt 



d k d k ® = p hW 2 -[p+^- )-{b t f (6) 



where W = -, 1 , D = pW, E = eW and Si = phW 2 Vi are the Lorentz boost factor, 

auxiliary variables correspond to density, energy, and momentum, respectively. Eq. ([2]), 
Eq. ([3]) and Eq. (jl]) represents the mass, energy, and momentum conservations. L u in 
the right hand side of Eq. ([3]) is a total neutrino cooling rate determined by microphysical 
processes which will be later explained. In Eq. (Tj0) it is noted that the relativistic enthalpy, 
h = (1 + e/p + p/p+ \b\ 2 j p) includes magnetic energy. Eq. ([5]) is the induction equation for 
the magnetic fields. B l are related to that in the rest frame of fluid as B % = Wb % — Wb t b l . 
Here b t is a time component of the 4-vector, 6 M . Eq. ([6]) is the Poisson equation for the 
gravitational potential, 

This newly developed code is an Eulerian code based on the finite-difference method. 
The numerical approach for solving the basic equations of ([2]), ([3]), and (jl]), consists of the two 
steps, namely, the transport and the source step. These procedures are essentially the same as 
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those of ZEUS-2D flStone fc Norman 



scheme of Van Leer is implemented ( 



1992 ). At the t ransport step, the second order upwind 



van Leerlll977l ). To handle the numerical oscillations, 



we employ an artificial viscosity. In the special relati vistic treatments, many forms for the 
compression heating are possible (IHawley et al.lll984bl ). In our code, we employ the form of 
■^Ldiy/^yWv 1 as the compression heating, which becomes the well-known artificial viscosity 
of von Neumann and Richtmyer under the Newtonian approximation. While not explicitly 
included in the above expression for the enthalpy, the contribution from the compression 
heating on the inertia is included in our calculations. The detailed status on the shock 
capturing using this term is shown at Appendix [Dl 

The time evolution of the magnetic fields is solved by induction equation, Eq. ([5]). 
In so doing, the code utilizes the so-called constrained transport method, which ensures 
the divergence free (V • B = 0) of the numerically evolved magnetic fields at all times. 
Furthermore, the method of characteristics (MOC) is implemented to propagate accurately 
all modes of MHD waves. The detailed explanation and the numerical tests are delivered in 
the appendix [B] The self-gravity is managed by solving the Poisson equation, Eq. (J6]) with 
the incomplete Cholesky decomposition conjugate gradient method. 

Together with these hydrodynamic procedures, the following microphysical processes are 
impleme nted in this code. We approximate the neutrino transpo rt by a multiflavor leakage 
scheme (jEpstein &: Pethickl Il98ll ; iRosswog fc Liebendorferl 120031 ) , in which three neutrino 
flavors: electron neutrino, v e \ electron antineutrino, v e \ and the heavy-lepton neutrinos, u^, 
Up, u T , v T (collectively referred to as ux), are taken into account. The neutrino reactions 
included are electron ca pture on proton and free nuclei; positron capture on neutron; ph oto-, 
pair, plasma processes ([Fuller et al.lll985l : iTakahashi et al.l[l978l : lltoh et al.lll989l . Em . We 



added a transport equation for the lepton fraction Yi(=Y e — Y e + + Y Ve — Y Pe ), 



dt 



-li 



(7) 



to treat their change due to the relevant charged current reactions, whose reaction rates 
are collectively represented by ji here, with Y e , Y e +, Y Ue , Yo e , 7/ being electron, p o sitron , 
electron neutrino, anti-electron n eutrino fraction, resp ectively (see Epstein fc Pethickl ( 1l98ll ); 
Rosswog fc Liebendorferl (120031 ); iKotake et al.l (120031 ) for details of the estimation of 7/). L u 
in Eq. (J3j) represents the total neutrino cooling rate which is also estimated by the scheme. 
As for the e quation of state ( EOS), we employ a realistic one based on the relativistic mean 
field theory (iShen et al.lll998l ). Since the pressure is not represented as the analytic function 
of density and internal energy like in the case of polytropic EOS, an iterative algorithm are 
employed to update the fundamental variables (see Appendix [C] for detail). 



In our two dimensional simulations, the spherical coordinate is used with 300 (r) x 60(0) 
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grid points to cover the computational domain. Axial symmetry and reflection symmetry 
across the equatorial plane are assumed. The radial grid is nonuniform, extending from to 
4.0 x 10 8 cm with finer grid near the center. The finest grid is set to 10 5 cm. The polar grid 
uniformly covers from 9 = to 9 = |. This choice of the grid numbers is sufficient for the 
aim of this paper as will be discussed in section 5. 

Finally we summ arize the difference on the numerical approach from our previous work 



(ITakiwaki et al.l 120041 ) . Most major development is the fully special relativistic treatment 
on magneto-hydrodynamics. And for the microphysical parts the cooling terms by neutrino 
contains contributions from not only u e but also z/ e and v%. These advances provide more 
reliable results on the magneto-rotational core-collapse. 



3. Initial Models 

We make precollapse models by taking the profiles of density, inte rnal energy, and 
electro n fraction distribution from a rotating presupernova model of E25 by lHeger fc Langer 
( 120001 ). This model has mass of 25M at the zero age main sequence, however loses the 
hydrogen envelope and becomes a Wolf Rayet star of 5.45 M before core-collapse. Our 
computational domain involves the whole iron-core of 1.69M . It is noted that this model 
seems to be a good candidate as a progenitor of the GRB since the lack of the line spectra 
of the ejected en yelopes ar e recon ciled with the observations of the supernovae associated 
with GRBs fe-g- lMeszarol tod )). 



Since little is known about the spatial dis tributions of the rotation and the magnetic 
fields in the evolved massive stars (see, however. ISpruitl (120021 )). we add the following rotation 
and magnetic field profiles in a parametric manner to the non-rotating core mentioned above. 
For the rotation profile, we assume a cylindrical rotation of 



n(x, z) = fi c 



7A 



7\ ' 



'x 2 + x$z* + z 4 

where Q is the angular velocity and X and Z denote distance from the rotational axis and 
the equatorial plane. We adopt values of the parameters, X and Z , as 10 7 cm, 10 8 cm, 
respectively. The parameter, X represents the degree of differential rotation . We assume 
the strong differential rotation as in our previous study (ITakiwaki et al.ll2004l ). 



As for the initial configuration of the magnetic fields, we assume that the field is nearly 
uniform and parallel to the rotational axis in the core and dipolar outside. For the purpose, 
we consider the following effective vector potential, 



A- = A 



0. 



(9) 
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Bn 



-r sin#, 



(10) 



2 r d + 

where A r ^^ is the vector potential in the r, 9, (ft direction, respectively, r is the radius, r is 
the radius of the core, and B is the model constant. In this study, we adopt the value of r 
as 2 x 10 8 cm which is approximately the size of the iron core at a precollapse stage. This 
vector potential can produce the uniform magnetic fields when r is small compared with 
r , and the dipole magnetic fields for vice versa. Since the outer boundary is superposed 
at r = 4 x 10 8 cm, the magnetic fields a re almost uniform in the computational domain as 
the previous work (ITakiwaki et al.l 12004 ) . It is noted that this is a far better way than the 
loop current method for constructing the dipole magnetic fields ( jSymbalistyl Il984l ) . because 
our method produces no divergence of the magnetic fields near the loop current. We set the 
outflow boundary conditions for the magnetic fields at the outer boundary of the calculated 
regions. 

We compute 9 models changing the total angular momentum and the strength of mag- 
netic fields by varying the value of Qq and B . The model parameters are shown in Table [TJ 
The models are named after this combination, with the first letters, B12, Bll, BIO, repre- 
senting strength of the initial magnetic field, the following letter, TW4.0, TW1.0, TWO. 25 
representing the initial T/IV4 7 !, respectively. Here T/|W| indicates the ratio of the rota- 
tional energy to the absolute value of the gravitational energy. The corresponding values 
of Q are 151rad/s, 76rad/s, 38rad/s for TW4.0, TW 1.0, TW0.25, respective ly. It is noted 
that the value of T/|VK| is 0.15% of the progenitor by lHeger fc Langerl (120001 ) and also that 
the specific angular momenta ranges from 0.5 to 1.5 jiQ for TWO. 25 to TW4.0 models with 
7ib = 10 16 cm 2 s -1 , which are in good agreement with the requirement of the collapsar model 
(IMacFadyen &: Woosleylll999l ). Current stellar evolution calculations predict that the rapidly 

er metalicity experien ces the so-called chemically homoge- 



rotating massive stars with smal 
neous cores during its evolution (jYoon &: Langer 



20051 ). Such stars are co nsidered to satisfy 



the r equirements of the collapsar model, namely rapid rot ation of the core (IWooslev fc Heger 
20061 ). According to a GRB progenitor model of 350B in lWoosley fc Hegerl (120061 ). the mag- 
netic field strength of the core reaches to ~ 10 12 G and the specific angular momentum is 
the order of j 16 ~ 1), with which our choices for the initial magnetic field and the initial 
rotation rate are reconciled. 
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Table 1: Models and Parameters 







T/\W\(%) 






0.25% 


1.0% 


4.0% 


10 10 G 


B10TW0.25 


B10TW1.0 


B10TW4.0 


B (Gauss) 10 n G 


B11TW0.25 


B10TW1.0 


B11TW4.0 


10 12 G 


B12TW0.25 


B10TW1.0 


B12TW4.0 



Note. — Model names are labeled by the initial strength of magnetic fields and rotation. T/|W| represents 
the ratio of the rotational energy to the absolute value of the gravitational energy. The corresponding values 
of fi in Eq. © are 151rad/s, 76rad/s, 38rad/s for TW4.0, TW1.0, TWO. 25, respectively. B represents the 
strength of the poloidal magnetic fields (see Eq. lflO)) ). The corresponding values of E m /\W\ is 2.5 x 10~ 8 , 
2.5 x 1CT 6 and 2.5 x 1CT 4 for 10 10 G, 10 11 G and 10 12 G, respectively with E m being the magnetic energy. 
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4. Results 

4.1. Hydrodynamics before Core-Bounce 

First of all, we briefly mention the dynamics before core bounce, when the gross features 
are rather similar among the computed models. The characteristic properties are summarized 
in Table H 

The story before core-bounce is almost the s a me as the canonical core-collapse super- 



novae with rapid rotation (see, e.g., iKotake et al.1 (120031 )). The core begins to collapse due 



to electron captures and the photodissociation of the iron nuclei, and eventually experiences 
the bounce at the subnuclear density by the additional support of the centrifugal forces. In 
fact, the central densities at bounce becomes smaller and the epoch till bounce is delayed as 
the initial rotation rates become larger (see pbnc and tbnc in Table [2]). 

As the compression proceeds, the rotational energy increases and reaches a few 10 52 erg 
at the moment of the bounce (seen from T/| W|bnc x WW in Table [2]). Given the same initial 
rotation rates, the values of T/|W|bnc do not depend on the initial field strength so much. 
This means that the angular momentum transfer is negligible be fore bounce, which is also the 



case of the Newtonian hydrodynamics (jYamada fc Sawail 12004 ). At bounce, the unshocked 



core becomes more flattened as the initial rotation rate becomes larger (compare panels in 
Figured]). The central protoneutron stars rotate very rapidly reaching to ~ 3000 rad/s 
with the typical surface magnetic fields of ~ 10 13 G to ~ 10 15 G for B10 and B12 models, 
respectively. From the table, it is also seen that the amplification rates of the magnetic fields 
(A amp ) are mainly determined by the initial rotational rates. One exception is the model 
B12TW4.0. Due to very rapid rotation with the highest magnetic fields initially imposed, 
the model bounces predominantly due to the magnetic force. As a result, the core bounce 
occurs earlier with the lower central density with less gravitational energy of the inner core 
than the models with the same initial rotation rate (see Table [2]). This earlier magnetically- 
supported bounce leads to the suppression of the amplification rate, which is exceptionally 
observed for this model. 

In this way, the hydrodynamic properties before bounce are mainly governed by the 
differences of the initial rotation rates. On the other hand, the differences of the magnetic 
field strength begin to play an important role on the dynamics later on. We will mention 
them in detail from the next sections. 
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Table 2: MHD properties till core bounce 



Model 


T/\W\ hQC 


Pbnc 


|Wbnc| 


E m /\W\i 


ni 


E m /\W\ hnc 


tip 
E m 


A 

^amp 


^bnc 


Names 




[10 14 g/cm 3 ] 


[10 53 erg] 












[ms] 


B12TW0.25 


0.10 


2.1 


1.1 


2.5 x 10" 


-4 


1.0 x 10~ 3 


0.3 


100 


245 


B11TW0.25 


0.10 


2.1 


1.1 


2.5 x 10" 


-6 


1.0 x 10~ 5 


0.3 


100 


245 


B10TW0.25 


0.10 


2.1 


1.1 


2.5 x 10" 


-8 


1.0 x 10~ 7 


0.3 


100 


245 


B12TW1.0 


0.18 


1.3 


1.1 


2.5 x 10" 


-4 


9.0 x 10~ 3 


0.07 


720 


295 


B11TW1.0 


0.18 


1.3 


1.1 


2.5 x 10" 


-6 


7.0 x 10~ 5 


0.07 


610 


295 


B10TW1.0 


0.18 


1.3 


1.1 


2.5 x 10" 


-8 


7.0 x 10~ 7 


0.07 


610 


295 


B12TW4.0 


0.20 


0.095 


0.68 


2.5 x 10" 


-4 


20 x 10~ 3 


0.3 


800 


477 


B11TW4.0 


0.19 


0.11 


0.74 


2.5 x 10" 


-6 


29 x 10~ 5 


0.1 


4400 


484 


B10TW4.0 


0.19 


0.11 


0.74 


2.5 x 10" 


-8 


31 x 10~ 7 


0.1 


4400 


484 



Note. — Characteristic properties before core bounce. T/|W|bnc is the rotational energy per gravitational 
energy at bounce, pbnc is the maximum density at bounce. E m /\W\ |i„; and -E m /|!^||bnc is the magnetic 
energy per the gravitational energy initially and at bounce, respectively, jf- is the ratio of the poloidal 
magnetic energy to the total magnetic energy at bounce. A amp represents the amplification rate of magnetic 

dof 

energy until core bounce, which is defined as A amp = (i? m |bnc)/(-E'm|ini)- ^bnc represents the time till bounce. 
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-3e*07 -lcH'CTT 1«*0 2**0 -2**07 -lr+07 Ic+0 2e+0 




-3e+07 -le+07 l*+0 2*+0 -2r+07 -li*07 l*+0 2c<0 



Fig. 1. — Snapshots for models B10TW0.25 (left) and B10TW4.0 (right), showing the degree 
of the rotational flattening at core bounce. In each panel, contour of density [g/cm 3 ] (left) 
and entropy per baryon [fc#] (right) is shown. Flattening of the unshocked core (the central 
low entropy region seen in the right side of each panel) is remarkable for the right panel. 
Note that the unit of the horizontal and the vertical axis is in cm. 
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4.2. Prompt vs. Delayed MHD Exploding Model 

After bounce, we can categorize the computed models into two groups, by the criterion 
whether the shock generated at bounce promptly reach the surface of the iron core or not. 
For later convenience, we call the former and the latter models as prompt and delayed MHD 
exploding model, respectively throughout the paper. The models and the corresponding 
groups are shown in Figure [2j To begin with, we choose typical model from the two groups 
and mention their properties in detail. 





T/\W\(%) 
0.25% 1.0% 4.0% 


10 10 G 

S (Gauss) 10 n G 
10 12 G 


122 ms 
72 ms 

32 ms 


96 ms 104 ms 


27 ms 32 ms 
20 ms 25 ms 



Fig. 2. — Classification of the computed models into the prompt (red blocks) or delayed 
(green blocks) MHD exploding model by the difference of tioookm shown in this table, which 
is the shock-arrival time to the radius of 1000 km after bounce. 



Prompt MHD Exploding Model The models classified into this group have strong 
magnetic fields and rapid rotation initially. Figure [3] shows the dynamics near core bounce. 
As seen from the bottom left panel, the shock at core bounce stalls in the direction to the 
equatorial plane at ~ 1.4 x 10 7 cm promptly (~ 3 ms) after bounce. However the shock in 
the direction of the rotational axis does not stall and becomes a collimated jet (see top right 
and bottom right). The wound-up magnetic fields are an important agent to explain these 
properties. 

The magnetic fields for the promptly MHD models are strong enough to power the jet 
already at the epoch of bounce. That is clearly shown in the top left panel, showing that 
the "plasma /?" = p/f-, being the ratio of the matter to the magnetic pressure, outside 
the unshocked core near the poles becomes very low (typically 10~ 2 ). From the right side 
of the bottom left panel, the toroidal magnetic field strength there reaches over 10 15 G. The 
dynamics around the poles are strongly affected by these strong magnetic fields. 

The three dimensional plots of Figure H] are useful to see how the field wrapping occurs. 
From the top left panel, it is seen that the field lines are strongly wound around the rotational 
axis. The white lines in the top right shows the streamlines of the matter. A fallback of 
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the matter just outside of the head of the jet downwards to the equator (like a cocoon) is 
seen. In this jet with a cocoon-like structure, the magnetic pressure is always dominant over 
the matter pressure (see the region where plasma j3 less than 1 in the right side of the top 
right panel of Figure [3D- This magneto-driven jet does not stall and penetrate to the surface 
of the iron core, which is essentially the repr o ducti on of the pioneering results in th e MHD 
supernova simulations by iLeBlanc fc Wilson! (119701 ) and its analysis by iMeier et al.l (119761 ). 
The speed of the head of the jet is mildly relativistic of ~ 0.3c, with c being the speed of 
light (the right side of bottom right panel of Figure [3]). At 20 ms after bounce, the jet finally 
reaches the surface of the iron core of ~ 10 8 cm. At this moment, the explosion energy, which 
will be a useful quantity for comparing the strength of the explosion among the models later, 
erg. 



reaches 1.4 x 10 50 
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Fig. 3. — Time evolution of various quantities characterizing the dynamics near after bounce 
for a prompt MHD exploding model (see text in section 4.2 of Prompt MHD exploding 
models). This is for model B12TW1.0. In each panel, the left side represents the logarithm 
of density [g/cm 3 ]. Time in each panel is measured from the epoch of bounce. At the top 

panels, the right side is the logarithm of the "plasma /3" = p/f-j indicated by "Beta". At 
the bottom left panel, the right side is the logarithm of toroidal component of the magnetic 
fields [G] , indicated by "-B p hi" • At the bottom right panel, the right side is the radial velocity 
in unit of the speed of light : c. Note that the unit of the horizontal and the vertical axis of 
all panels are in cm. 
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Fig. 4. — Three dimensional plots of entropy with the magnetic field lines (left) and the 
streamlines of the matter (right) during the jet propagation for models of B12TW1.0 (top) 
and B10TW1.0 (bottom), at 20 ms and 94 ms after bounce, respectively. The outer edge of 
the sphere colored by blue represents the radius of 7.5 x 10 7 cm. Note that the model of the top 
and bottom panel belongs to the prompt and delayed MHD exploding model, respectively. 
These panels highlight not only the wound up magnetic field around the rotational axis 
(left), but also the fallback of the matter from the head of the jet downwards to the equator, 
making a cocoon-like structure behind the jet (right). 
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Delayed MHD Exploding Model The models with weaker initial magnetic fields belong 
to the delayed MHD exploding model (see Figure [2]). In the following, we explain their 
properties taking model B10TW1.0 as an example. It is noted again that this model has the 
same initial rotation rate with model B12TW1.0 of the previous section, but with the two 
orders-of-magnitudes weaker initial magnetic fields. 

In the case of model B10TW1.0, the shock wave at bounce stalls in all directions at 
~ 1.5 x 10 7 cm. As shown in the top left panel of Figure the plasma (3 is so high that 
the magnetic fields play no important role before bounce. After the shock stalls, the stalled 
shock begins to oscillate. The middle left and the bottom left panel shows the prolate and 
oblate phase during the oscillations, respectively. Until ~ 70 ms after bounce, the oscillation 
of the shock front continues diminishing its amplitude. Approximately the number of the 
oscillations is about 5 times this time. Without the magnetic fields, the oscillation should 
cease settling into the equilibrium state with the constant accretion through the stalled-shock 
to the center. However during this oscillation, the magnetic fields behind the stalled shock 
gradually grow due to the field wrapping and the plasma (3 around the polar regions becomes 
low as seen from the right side of the top right panel. Soon after the toroidal magnetic fields 
become as high as ~ 10 15 G behind the stalled shock (see the middle right panel), the stalled 
shock near the pole suddenly begins to propagate along the rotational axis and turns to be a 
collimated jet (see the bottom right panel). This revived jet does not stall in the iron-core. 
This is the reason why we call this model as the delayed MHD exploding model. The speed 
of the jet reaches about 5.5 x 10 9 cm/s (see the bottom right panel). Also in this jet, the 
toroidal component of the magnetic fields is dominant over the poloidal one and a fallback 
of the matter is found in the outer region of the jet (cocoon) as in the case of the promptly 
MHD exploding model (see the bottom two panels of Figure H]). At ~ 96 ms after bounce, 
the jet reaches ~ 10 8 cm. The explosion energy at that time reaches 0.094 x 10 50 ergs. 

As mentioned, the dynamical behaviors between the prompt and delayed MHD ex- 
ploding models after bounce seem apparently different. However there are some important 
similarities between them, which we discuss from the next section. 
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Fig. 5. — Same as Figure 3 but for the quantities showing the dynamics near after bounce for 
a delayed exploding model (see text in section 4.2:Delayed MHD exploding model). This is 
for model B10TW1.0. During the oscillations of the stalled shock after bounce (from middle 
left to bottom left), the magnetic fields behind the stalled shocks become large enough, due 
to the field wrapping (top and middle right), leading to the shock- revival for the formation 
of the magnetically-dominated jet (bottom right). Note that the unit of the horizontal and 
the vertical axis of all panels is in cm. 
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4.3. Similarities of Prompt and Delayed MHD Exploding Model 

In this section we focus on the similarities between the prompt and delayed MHD 
exploding models. 

From Figure it can be seen that the radial velocities and the magnetic fields of the jets 
are quite similar among the models regardless of the prompt or delayed exploding models. 
Typical values of the toroidal magnetic fields are 10 14 — 10 16 G and typical velocities are 
10 — 30% of the speed of light. The opening angles of the jets are also similar. The width 
of this jet is about 8 x 10 6 cm when the jet reaches 7.5 x 10 7 cm, which means that the half 
opening angle of the jets is about 6° at this time. These characteristic values of the jets are 
summarized in Table [3j 

Detailed properties of the jets in the vicinity along the rotational axis are shown in 
Figures [7J and [H] to see the origin of these similarities. We fix the initial rotation rate and 
the initial field strength in Figure [7J and [HI respectively, to see their effects separately. In 
Figure UJ, the initial rotation rate is T/|VT| = 1.0% and the different lines correspond to the 
difference between the initial magnetic fields from 10 12 (B12) to 10 10 G (BIO). In Figure [BJ 
the initial magnetic field is 10 11 G and the different lines corresponds the difference in the 
initial rotation rates. 

From the top and middle panels of Figures [7J and [HI we find that the radial profiles 
of the toroidal magnetic field, the plasma (3 (0.1 — 0.01), the density, and the velocity, 
are rather similar behind the shock whose position can be seen from the discontinuity at 
~ 700 km. Above all, it is surprising to see the remarkable similarity in the profiles of 
the toroidal magnetic fields behind the shock among the models (top left in Figures [7J and 
[H]). The typical strength behind the shock is seen to be ~ 10 15 G. This critical strength of 
the toroidal magnetic field for the shock-revival is estimated as follows. The matter behind 
the stalled-shock is pushed inwards by the ram pressure of the accreting matter. This ram 
pressure is estimated as, 

P = 4 x 10 28 ( J, \ ( ^— rV erg/cm 3 , (11) 

Vio lo g/cm 3 y V 2 x 10 9 cm/ S y s/ v ; 

where the typical density and the radial velocity are taken from Figure [1] and the bottom 
right panel of Figure [5], respectively. When the toroidal magnetic fields are amplified as large 
as ~ 10 15 G due to the field wrapping behind the shock, the resulting magnetic pressure, 
can overwhelm the ram pressure, leading to the magnetic shock-revival. The origin of the 
similarity of the jets seen in Figure [3] comes from this mechanism. We find that this process 
works in all the computed models. It is noted that the importance of the magnetic -shock 
revival was noticed also in the analytic models by lUzdensky Sz MacFadyen j2007cM . In 
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addition to their expectations, our simulations show that the explosion energy becomes 
smaller than their estimations because the magnetic tower cannot be wider as they assumed. 

From the bottom panels of Figure [HI it can be seen that the poloidal fields behind the 
shock front do not depend on the initial rotation rate so much given the same initial field 
strength, while the difference of the poloidal magnetic fields behind the shock in the bottom 
panels of Figure [7] simply comes from the difference in the initial field strength. This feature 
is regardless of the prompt or delayed models. 
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Fig. 6. — Contour of radial velocity (left side) and toroidal magnetic fields (right side in 
each panel) for various models, showing the similarity of the jets between the prompt and 
delayed MHD exploding models. The right and left panels correspond to the rapidly rotating 
(T/|VK| = 4.0%) and the slowly rotating T/|VK| = 0.25% models, respectively. From top to 
bottom panels, the initial strength of the magnetic fields changes from strong (B12 models) 
to weak (BIO models). The time from bounce is shown in the top right part of each panel, 
indicating the difference between the prompt and the delayed models. Note that the unit of 
the horizontal and the vertical axis of all panels is in cm. 
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Fig. 7. — Various profiles of jets along the rotational axis, for models B12TW1.0 at 20 ms 
(solid), B11TW1.0 at 24 ms (short-dashed) and B10TW1.0 (dotted) at 94 ms after bounce, 
respectively, showing the effects of the difference of the initial magnetic fields (fixing the 
initial rotation rate of T/|W| = 1.0%). Top left and right is the toroidal magnetic field (B^) 
and plasma beta {(3 = P/j-). Middle left and right is the radial velocity [cm/s] and the 
density [g/cm 3 ]. Bottom left and right is the poloidal magnetic field (B p ) and the ratio of 
toroidal to poloidal magnetic field {B^/B p ). Note that the shock position is approximately 
700 km as seen from the discontinuity of these profiles. It should be noted that the toroidal 
fields here are not just at 9 = 0° (the fields are zero there) but at the closest mesh to the 
axis of 6 = 1.5°. 
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Fig. 8.— Same as Figured but for models B11TW4.0 at 28 ms (solid), B11TW1.0 at 24 ms 
(short-dashed), and B11TW0.25 (dotted) at 64 ms after bounce, respectively, showing the 
effects of the difference of the initial rotation rates (fixing the initial field strength of 511). 
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Table 3: Characteristic properties of jets 



Name 


^-exp 


^lOOOkm 




E 

-^ CX P 1000km 




[ms] 


[ms] 


[cm/s] 


[10 50 ergs] 


B12TW0.25 





32 


6.0 x 10 9 


1.3 


B11TW0.25 


48 


72 


3.7 x 10 9 


0.05 


B10TW0.25 


92 


122 


3.8 x 10 9 


0.02 


B12TW1.0 





20 


8.0 x 10 9 


1.4 


B11TW1.0 


10 


27 


8.0 x 10 9 


0.23 


B10TW1.0 


70 


96 


5.5 x 10 9 


0.094 


B12TW4.0 





25 


6.0 x 10 9 


1.0 


B11TW4.0 


16 


32 


6.0 x 10 9 


0.10 


B10TW4.0 


65 


104 


5.5 x 10 9 


0.006 



Note. — Properties of the jets. i cxp is the duration from bounce to the revival of the stalled shock 
due to the field wrapping, tioookm represents the duration required for the jet to reach ~ 1000km after 
bounce. Explosion energy, ^expjoookm' an< ^ ^he velocity; Ujet is estimated at the moment. ^expioookm are 
normalized as 10 50 erg. For the definition of the explosion energy, see Eq. (|12p. 
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4.4. Dependence of Jet Arrival Times and Explosion Energies on Initial 
Rotation Rates and Magnetic Field Strengths 

In the previous section, we discuss the similarity among the computed models. From 
this section, we move on to discuss the differences among them. 

Jet Arrival Time First we discuss the "jet arrival time" shown in Table [2J which is the 
timescale when the jet reaches the outer edge of the iron core of ~ 1000km. As discussed in 
the previous section, this timescale is mainly determined how long it takes for the magnetic 
fields behind the shock to become as large as the critical toroidal magnetic fields (~ 10 15 G) 
SIS db result of the field wrapping. 

From the top left panel of Figure [H it is seen that the strong initial magnetic fields 
shorten the jet arrival time. This tendency is seen in all the computed models regardless 
of the prompt or delayed exploding models. When the initial magnetic fields are strong 
enough (~ 10~ 4 of the gravitational energy), the jet arrival times between the different 
initial rotational models become almost the same. In this case, the critical magnetic fields 
for the shock-revival are already generated by the compression before core bounce. So the 
strong magneto-driven jets can produce the prompt MHD explosions in a similar way. For 
the rapidly rotating models (the sequence of TW1.0 and TW4.0), it is seen that the decrease 
in the rate of the jet arrival time as a function of the initial E m /\W\ becomes smaller when 
the initial E m /\W\ is larger than ~ 10 -6 (see the kink in the panel). This is because too 
strong magnetic fields transport the angular momentum of the protoneutron star outwards, 
leading to the suppression of the efficiency of the field wrapping after bounce. 

In the top right panel of Figure [HI the dependence of the jet arrival time on the initial 
rotation rate is shown. By intuition, the jet arrival time may become shorter as the initial 
rotation rates become larger since the field- wrapping should become more efficient. The 
panel shows that this is true for moderately rotating models of the initial T/|W| less than 
0.01, but not true for the more rapidly rotating models. This can be explained as follows. 
Too rapid rotation of the core hinders the central core from collapsing due to the stronger 
centrifugal forces. This feature is clearly shown in the middle left panel of Figure [9] showing 
the density profiles. The density near the center is ~ 100 times lower than that for the slowly 
rotating models. Since the angular momentum is well conserved before bounce (see section 
4.1), the inner core (< 20 km) gains smaller angular velocities for rapidly rotating models 
by the weakened compression as seen in the middle right panel of Figure M Reflecting these 
aspects, the amplification rate of the magnetic fields (^j^/Em) near core-bounce becomes 
smaller for the most rapidly rotating model (TW4.0) as seen from the bottom panel of Figure 
9. This suppression makes the jet arrival time almost constant or longer as the initial T/|VF| 
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becomes larger than ~ 0.01 as in the top right panel of Figure [9j 

Explosion Energies In addition to a wide variety of the jet arrival times, we find a large 
difference in the strengths of the magnetic explosions. 

As a measure of the strength, we define the explosion energy as, 



here ei oca i is the sum of e^ in , e int , e mag and e grav , with being the kinetic, internal, magnetic, and 
gravitational energy, respectively (see Appendix I A . 1 1 for their definitions in special relativity) 
and D represents the domain where the local energy is positive, indicating that the matter 
is not bound by the gravity. The explosion energy is evaluated when the jet arrives at the 
radius of 1000km at the polar direction. The value of the explosion energy is summarized 
in Figure [10J Generally speaking, it is found that the explosion energies becomes larger for 
the prompt MHD exploding models (red) than the delayed MHD exploding models (green). 

What makes the difference on the explosion energies at the shock breakout from the iron 
cores? Firstly, the initial strength of the magnetic field is the primary agent to affect the 
explosion energies. The explosion energies are larger for models with the larger initial fields 
as seen in Figure 11. Secondly, the geometry of the jets has also effects on the explosion 
energies. Figure [11] shows the toroidal magnetic fields (left side) and the local energy (right 
side) in the jets from the stronger to the weak magnetic fields models (from top to bottom 
panels) at the shock breakout. In each right panel, it is noted that the regions with the 
positive local energies (ei oca i > in Eq. (12)) are drawn with color scales and the regions 
with black are for the regions with the negative local energies. It is seen that the regions 
where the local energy is positive mostly coincide with the regions where the strong toroidal 
magnetic fields are generated. As the initial field strength becomes larger, the regions where 
the local energy becomes positive, becomes larger (i.e. the jets become wider), leading to 
the larger explosion energies. In the case of the delayed exploding model (the right panel 
in Figure [TTj) . it is found that the width of the jets becomes narrower, which results in the 
smaller explosion energies. Although the properties of the jets just on the rotational axis 
are similar among the models seen from Figures 7 and 8, the lateral structures of the jets is 
found to have influence over the explosion energies. 
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Fig. 9. — Dependence of the jet arrival time on the initial magnetic field strength (top left) 
and the initial rotation rate (top right). Here E m /\W\ and T/|W| represent the ratio of 
initial magnetic and rotational energies to the gravitational energy, respectively. The jet 
arrival time is the duration for the jet to reach the outer edge of iron core of ~ 1000km. 
Middle left and right panel shows the distribution of the density and the angular velocity for 
models B11TW4.0, B11TW1.0 and B11TW0.25 at core bounce as a function of the equatorial 
radius, respectively. Bottom panel shows the temporal amplification rate (^^-/E m ) near core 
bounce. 
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Fig. 10. — The explosion energy for all the models. The energies are normalized to 10 erg. 
Note again that the red and green blocks indicate the prompt and delayed MHD exploding 
models, respectively. See Eq. ( fT2l) for the definition of the explosion energy. 
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Fig. 11. — Contour of toroidal magnetic field and the local energy for models B12TW1.0, 
B11TW1.0 and B10TW1.0 near the shock breakout from the iron core. In each panels, the 
logarithm of toroidal magnetic field [G] (left side) and the logarithm of local energy [ergs/cm 3 ] 
(right side) is shown, respectively. It is noted for each right panel that the regions with the 
positive local energies (ei oca i > in Eq. (12)) are drawn with color scales and the regions with 
black are for the regions with the negative local energies. The timescales given is measured 
from the core-bounce indicate the prompt or delayed MHD exploding models. The unit of 
the horizontal and the vertical axis of all panels are in cm. 
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5. Summary and Discussion 



We performed a series of two-dimensional MHD simulations of rotational core-collapse 
of magnetized massive stars. The main motivation was to clarify how the strong magnetic 
fields and the rapid rotation of the core affect the magnetic explosions. To handle the very 
strong magnetic fields, we developed a new code under the framework of special relativity. A 
novel point is that the microphysics such as the realistic equation of state and the neutrino 
cooling are implemented to the special relativistic MHD code. Due to these advantages, 
our computation can achieve a longer time-stretch of the evolutions compared to previous 
studies. The obtained results can be summarized as follows. 



• Magnetically powered jets are commonly found in all the computed models. In the 
jets, the magnetic fields are dominated by the toroidal components as a result of 
the field wrapping. For the profiles and strengths of the toroidal fields behind the 
jets, we find a remarkable similarity. We find that the jet-like explosions occur when 
the magnetic pressure behind the shock becomes strong, due to the field wrapping, 
enough to overwhelm the ram pressure of the accreting matter. The required toroidal 
magnetic fields are similar of ~ 10 15 G, which can be also understood by a simple order- 
of-magnetite estimation. Reflecting the similarity in the mechanism of producing jets, 
global properties of the jets such as the velocities (~ 20% of the speed of light) and the 
half opening angle of the jets (~ 6°) are also found to be similar among the computed 
models. 

• The timescale before the onset of the magnetic shock-revival are quite different depend- 
ing on the initial strengths of rotation and magnetic fields. When the initial strengths 
of rotation and magnetic field are larger, the jet can be launched just after the core- 
bounce, which we called as the promptly MHD exploding models. We furthermore 
find that even for the model with the weaker initial field and slow rotation, the jet-like 
explosions can occur after sufficient field wrapping to reach the critical field strength, 
which we called as the delayed MHD exploding models. In this case, the explosion 
can be delayed about ~ 100ms after bounce. The explosion energy also strongly de- 
pends on the time difference before the shock-revival. The stronger initial magnetic 
fields make wider exploding regions, leading to the larger explosion energy. The largest 
MHD-driven explosion energy obtained is ~ 10 50 erg. 



In addition to the magnetic shock-revival, the neutrino-driven shock revival, namely 
neutrino heating from the new ly born protoneutron star may energetize the jets as sug- 
gested by iMetzger et ajj (120071 ). Although we treated only the neutrino coolings in our 
computations, we try to estimate the effect of the neutrino heating in the following way 
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to see which one could be more important to produce the jets. We compare the energy 
gained by the neutrino heating to the magnetic energy. The specific neutrino heating rate 
due to neutrino absorptions [y e + n —>■ e~ + p and v e + p — > e + + n), which are the dom- 
inant heati ng processes at ce r tain r adius in the postbounce phase, can be estimated using 



Eq. (10) of iQian fc Wooslevl (11990 1. For model B11TW1.0 at 22 ms after bounce, the jet 
reaches ~ 5 x 10 7 cm and the density at the head of the jet, pj et , is ~ 1.5 x 10 8 g/cm 3 . At 
that time, the neutrino sphere locates at ~ 7 x 10 6 cm. The average energies and the lu- 
minosity for the electron and anti-electron neutrinos there are about 10 and 14 MeV, and 
50 and ~ 7 x 10 51 erg/s. In this setup, the heating rate due to the neutrino absorptions, 
q U N, reads ~ 7.7 x 10 22 MeV/s/g. In the same way, the heating rate due to the neutrino 
pair-annihilation (y + v —>■ e~ +e + ) can be estimated as, ~ 1.0 x 10 16 MeV/s/g, which is neg- 
l igible compared to the neutrino absorptions, albeit wi t h the general relativistic corrections 



( ISalmonson fc Wilson! Il999l ; lAsano fc Fukuyamal l200ll . |2000j), for far outside the neutrino 



spheres where the magnetic shock revival occurs. If the neutrino luminosity maintains dur- 
ing the characteristic time scale for the delayed MHD explosion, namely Allayed ~ 100 ms, 
the fluid element may gain q uN (pj et /1.5 x 10 8 g/cm 3 )(At delayed /100ms) ~ 1.6 x 10 25 erg/cm 3 
from the neutrino heatings. This estimation shows that the energy deposition from neutrinos 
is quite smaller than the magnetic energy behind the jet, ~ G ^ ~ 5.4 x 10 28 erg/cm 3 . 

Here, we omitted the neutrino cooling to maximize the effect of the neutrino heating. This 
simple estimation shows that for the models computed in this paper, the neutrino heatings 
could be minor compared to the magnetic effects. However for more weakly magnetized and 
slowly rotating cores, the neutrino heating may overwhelm the magnetic shock-revival. To 
seek the criteria which mechanisms can be dominant, is important, however, beyond scope 
of this study. 

With respect to our numerical computations, we have to give some discussions. First 
of all, we discuss the convergence of our numerical results. As mentioned in section [21 we 
have taken the grid numbers of N r = 300 and Ng = 60 as a fiducial value with N r and Ng 
being the grid numbers in radial and polar direction, respectively. For the convergence tests, 
we vary N r from 200, 400, to 1000 fixing N e = 60 and N e from 30, 90, to 120, while fixing 
N r = 300. Taking model B11TW1.0, we follow the dynamics for models with the different 
numerical resolutions till ~ 30 ms after bounce when the magneto-driven jets come out of the 
central iron cores. We pay attention to the (spatially-integrated) magnetic energies among 
the models, because they are a good indicator to see the degree of the amplification of the 
magnetic fields. Varying the coarser to the finer grid resolutions for the radial and polar 
direction, the relative change in the magnetic energies from the fiducial model is found to 
be in the range of +0.15% — 0.35% and ±3.0%, respectively. These resolution tests seem 
to support the convergence of the numerical values obtained in this paper. However we 
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shall note that better angular resolutions are needed to resolve the collimation of the jets 
further than ~ 4 x 10 7 cm from the center. This is because the intervals of the radial 
mesh, Ar, become logarithmically large in the spherical coordinates, which may lead to 
the overestimation of the collimation. Since the jets propagate rather cylindrically around 
the polar axis later on, we think it a good way to switch from the spherical coordinates to 
the cylindrical coordinates at certain radius to sustain better angular resolutions, however, 
beyond scope of this study. Next we estimate the numerical dissipation of the magnetic 
fields. During the infall epoch, the toroidal magnetic flux is a good quantity to see the flux 
conservation. The difference between the initial toroidal magnetic flux, <&- m i, and the flux at 
the bounce, $bnc ; is estimated to be A$ = $i "'^ #bnc ~ 1%. This means that the numerical 
magnetic dissipation is treated to be small for our numerical code. 

Then we move on to address a few imperfections in our simulations. Magneto-rotational 
Instability (MRI) has been pointed out to be active in the outer layer of protoneutron star 



which rotates strongly differentially (lAkiyama et al.ll2003l ; iKotake et al.ll2004al ; iMasada et al. 
20061 ). The wavelength of maximum growth rate of the linear instability, A = ^f 2 -, becomes 

1 /2 

10 Hq) ( 10 ^ cm ) km f Balbus fc Hawleylll998 ). where our numerical grid is 



5 (»g=i) (- 



10 10 g/cm 3 
P 

3km there. Since 10 — 100 times finer mesh than the wavelength is required for resolving 



the MRI (IShibata et al.l 120061 ) . our simulations are insufficient to take into account the field 
amplification due to the MRI. This is a very important task remained to be investigated, 
although the computational costs are still expensive. As for the microphysi c s, the neutrino 
heating is not included as mentioned above (see, however, iBurrows et al.l (120071 )). Since 
the timescale for the magnetic explosions in the polar direction are much earlier than the 
neutrino heating or the g-mode excited explosions, we think that both of them basically 
should play a minor role here. On the other hand, the neutrino heating could be helpful 
for producing the explosion in the direction of the equatorial plane where the field-wrapping 
induced magnetic explosions are unlikely to occur. Finally the general relativistic (GR) 
effects are not incorporated in our simulations. During our simulation time, the central 
protoneutron stars do not collapse to the black holes as inferred from a simple estimation 
of the Schwarzshild radius. Thus we think that the GR effects may not drastically change 
our results qualitatively. After the jet break-out from the star leaving behind the narrow 
funnel, the mass accretion in the direction of the equatorial plane may lead to the black 
hole formation, which provides us with the initial condition of the collapsar. This would be 
another interesting subject as an application of this study. 

Bearing these caveats in mind, we state some speculations based on our results. The 
protoneutron stars obtained here are with the poloidal magnetic fields of order 10 15 [G] and 
with its rotation period of an order of milliseconds, which could be the origin of the mag- 
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netar according to a hypothesis by iDuncan fc Thompson! (Il992l ). If so, it means that the 
magnetically-driven jets could be associated with the formation of the magnetars. Our results 
suggest that the toroidal component of the magnetic field is dominant in the young magne- 
tars. The large magnetic energy of the toroidal component stored in deep crusts and cores 
of the magnetars would be transported outside via Alfven waves and be released as giant 



flares of the SGRs (IDuncan fc Thompsonlll992l ; iThompson fc Duncanll200ll ; iThompson et al. 



20021 ). Some observational evidences support this picture. For example, quasi periodic oscil- 



from the Alfven wave from the interior of the star (IRea et al.ll2006l ; llsrael et al 



Israel et al. 


2005). There 


the XRF ( 


Mazzali et al. 



20061 : iMaeda et al]l2007al fbh . While the ordinary G RBs require the highly relativistic eiecta 



the mildly relativistic ejecta is favorable for XRFs (jToma et al. 2007 ; Soderberg et al. 



2006 



Ghisellini et al.l 120071 ). which may be the case here because the magneto-driven jets become 
only mildly relativistic due to the high baryon loading of the matter along the rotational 



axis. 



There remain more rooms to be investigated applying our simulations. In this study, we 
employed one progenitor model. Since the ac cretion rate of the matter to the protoneu tron 
stars should depend on the progenitor mass (IHeger et al.ll2005l ; IWoosley Hegerl 120061 ). we 
think it very important to investigate how the criterion of the magnetic shock-revival changes 
with the progenitor models. Moreover the initial configurations of the ma gnetic fields, which 



are st ill highly uncertain, could be changed in a systematic manner like in lSawai et al.l (12005 



20071 ). to see their effects on dynamics. While this study focused on the shock-propagation in 
the iron cores, in which the jets become only mildly relativistic, we plan to continue to follow 
the dynamics later on till the jets break out of the stars (phase 4 i n the introduction), in 
which t he jets are expected to b e relativistic. Very recent studies by iKomissarov fc Barkov 
( 120071 ); iBucciantini et al.l (120071 ) are on this line. Our simulation can be more consistent 
than their studies in the sense that we start the simulations from the onset of the core- 
collapse, and that the protoneutron stars are not excised like their models. By continuing 
the simulations of the jet propagations till the shock-breakout, we plan to study the possible 
connection between the magnetically-driven jets obtained here and the origins of the XRFs 
in the forthcoming work (Takiwaki et al. in preparation). 
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A. Derivation of the Basic Equations for SRMHD 



In this Appendix, we summarize formalisms on the basic equations and the numerical 
tests for our new l y dev el oped SRMHD code. For the formalisms, we follow the derivation of 
De Villiers etaJ fl2003h : lHawlev et al.l Jl984al lbT). For convenience, we proceed the derivation 
keeping the metric general forms, i.e., ds 2 = — a 2 dt 2 + 7;j (dx 1 + /3'dt) (dx^ + /3 J dtJ . where 
a is the lapse function, f3 is the shift vector and is the spatial 3-metric. And we take the 
Minkowski metric later. 



This preprint was prepared with the AAS IATgX macros v5.2. 
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There are four fundamental magnetohydrodynamic equations. The conservation of 
baryon number is 

dupU" = (Al) 

where p,£7 M (/i = 0, 1, 2, 3) are baryon mass density and 4- velocity at each point. The conser- 
vation of the stress-energy is 

= (A2) 

where T^ u is stress-energy tensor and Maxwell's equations 

dpF 1 ™ = 4vrJ", (A3) 
q* F iw = o. (A4) 

where F^ v is the antisymmetric electro-magnetic tensor and 

= -e^ Sa F Scr (A5) 



is dual of F^ v . Maxwell's equations are supplemented by the equation of the charge conser- 
vation d^jv = 0. 

The energy momentum tensor consists of perfect fluid parts and electromagnetic parts, i.e. 



= ph*U»U p + V f v + (f^F™ - ~F afi F af, g> H 'J 



(A6) 



where h* = (1 + e/p + p/p) is the relativistic enthalpy with e and p being the internal energy 
and the pressure, respectively. 

For later convenience, we define magnetic induction in the rest frame of the fluid, 

6" = -^=*F^U U . (A7) 

V47T 



We adopt the ideal MHD limit and assume infinite conductivity (the flux- freezing 
condition), where in the electric field in the fluid rest frame is zero, i.e., F^ U U U = 0. 

Combining Eq. (1A5|) with (Eq. IA7I) and conditions for infinite conductivity, we obtain 

= e a ^ v ^b a U p . (A8) 

The orthogonality condition 

= (A9) 
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follows directly from Eq. (IA7I) . 

The induction Eq. ()A4j) can also be rewritten by substituting the definitions, 

d a (U a b 13 - b a U p ) = 0. (A10) 

By expanding this equation using the product rule and applying the orthogonality condition 
Eq. (1A9[) . we obtain the identity 

U v 6 M V M U v = 0. (All) 
It is useful to rewrite the energy momentum tensor as 

T» v = (ph* + \b\ 2 )U tM U u + (p + ^-)g^ -Vb v . (A12) 

We have to expand basic equations in terms of the code variable, and transform the 
equation for auxiliary density, energy and momentum functions D = pW, E = eW, Si = 
phW 2 Vi. Finally the set of variables D, E, Si, Bi will be evolved through the basic equations 
transformed here. 

The equation of baryon conservation (lAlj) can be expanded in terms of the code variables 
easily, 

d t D + ^-d j (D y /yV i ) = 0. (A13) 
The equation of energy conservation is derived by contracting (1A2j) with U v , 
U V V^ V = U u vA (ph*+\\bf) U»U U + \P+^Y-j f v ' -Wb v \ =0. (A14) 



By using the identity (1A11I) and (lAljl . we obtain the local energy conservation 

V M (pe^) + PV M ^ = 0, (A15) 

Applying the definition for the auxiliary energy function E, the energy equation is rewritten 
as follows: 

d t (E) + ±- 9i ( ^7 E V 1 ) + P d t (W) + ^- di (V7 W V 1 ) = 0. (A16) 



The momentum conservation equations follow from 

V M T» u = vA(ph*+ \\b\\ 2 ) u»U v +[P + S-]5 li v -V i b v \=0. (A17) 
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This equation can be rewritten as 

.VI* 4- 

2a S* 




\bf\ I I 



^ a v^6^ 6, - - b a bp d v g a/3 = 0. (A18) 



" v ' 2 



To obtain the final form of the equations, multiply (IA18I) by the lapse a, split the p, index 
into its space {%) and time (t) components, and restrict v to the spatial indices (J) only: 

d t (Sj - a b b l ) +J-d iy fy (Sj V*-a b 3 U) 
y7 



\ {^ + ab ^) djg^-adJp+^fY (A19) 



The v index can be restricted to the spatial indices because the equation that arises from v = 
t for the time components of momentum and magnetic fields is redundant, corresponding to 
the total energy conservation equation. In our formalism, we solve the Eq. (1A16j) separately 
for the internal energy Taking the following metric, 

ft. = ( - 1 ^ ) • (A20) 

where 7^ is the spatial metric whose concrete description depends on coordinate system. 
Finally we describe our treatment of the gravity. Under the weak field limit, time-time 
component of the metric, q t t, takes the form of — (1 — 2$) where $ is Newtonian gravitational 



potential (e.g. IShapiro &: Teukolskylll983l ). The third term of the momentum equation then 
becomes, 

1 S a S 13 

- -—d j9a ^ phW'd^. (A21) 

Under this limit, Einstein equation becomes the Poisson equation for the gravitational po- 
tential (see Eq. (jHJ)). Since the origin of the source term is the tt component of the energy 
momentum tensor, we replace p in the ordinary Newtonian limit with T u - as we discussed 
in the appendix. 

The validity of using the Newtonian potential in the core-collapse simulations here may 
be discussed by seeing the value of compactness parameter, 2G ^ r ^ , where r is the radius and 
M(r) is the enclosed mass within r. In the vicinity of the protoneutron star of ~ 1.2M with 
the typical size of ~ 20km, the parameter is ~ 0.18. So the error caused by neglecting higher 
order metric perturbations is estimated to below ~ 3%. The qualitative features found in 
this paper may be unchanged due to the incursion of the GR. 
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A.l. Energy Descriptions 

We need to modify the description of energy from the Newtonian one to the special 
relativistic one. The total local energy, ei oca i, is defined by sum of the various energy: 

^local ^kin ^int ^mag ^grav (A22) 

where eki n , ei nt , e grav and e mag is kinetic energy, internal energy, gravitational energy and 
magnetic energy, respectively. Their specific description is as follows: 

e kin = P W(W-1), (A23) 

e in t = eW 2 +p(W 2 -l), (A24) 

e grav = -phW 2 $, (A25) 

6° 2 



e - = B \ 1 -m)-w^- (A26) 

These descriptions are used for the calculations of the explosion energy in Subsection 14.41 



B. Special Relativistic MOC 



The method of characteristics (MOC) is popularly used in the magneto-hydrodynamical 
simulations. In this algorithm the magnetic fields are evolved along the characteristic lines 
of the Alfven waves. Deta iled procedure for this algorithm for the Newtonian case is given 
Stone Norman! (Il992l ). For the special relativistic (SR) computations, we derive the 



m 



solutions of the SR Alfven waves in an analytic form, 



Dt 



Dt 



D 

Dt 

D 
Dt 



def 



dcf 



d 

dt + 







d 



(1 - 6V 

Va + ■ 



+ 



phW) dxj 
d 



phW 



dt (l + b t / y /p~hW)dx 



(B2) 



(B3) 



(B4) 



where W, p, and h is the Lorentz factor, density and enthalpy respectively. Vj and bj is the 
perpendicular component of the velocity and the magnetic field to the Xj directions. 



In the subroutine for solving SR MOC in the code, the velocity and the magnetic fields 
are updated at half-time step along the characteristics using the above equations. By giving 
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the analytic forms, it is readily seen that the speed of the propagation is guaranteed to be 
below the speed of light even for the regions where the density becomes low and the magnetic 
fields become strong, which is quite important for keeping the stable numerical calculations 
in good accuracy 



Alfven Wave Propagation The propagation of a liner Alfven wave is a basic test prob- 
lem of MHD simulation. We consider a constant background magnetic field, B x , and fluid 



velocity, v x . And we add small transverse perturbations with velocity, v z 



J y)i 



and magnetic 



field, B z (By). In this situation b l = J2v k bk 



Alfven wave becomes 



DWv z + b z /^ph 



Dt 

BWv z - b z j^ph 



Dt 



D 

Dt 

D 

Dt 



d v 



therefore the analytic solution for the 

(B5) 
(B6) 

(B7) 



'DhW 







dt 



V x + 



'DhW 



'DhW 



dx ' 




1 + v- 



'DhW 



dx 



If we take 



Wv z + b z /^ph = 



(B9) 



the minus mode does not propagate, we assume B x = 0.09, v x = 0.08 for the Newto- 
nian Alfven wave and B x = 0.9, v x = 0.8 for the relativistic Alfven wave. We take v z as 
10~ 7 v x sin(27rx) and B z is determined from b z in Eq. flB9[) . The result is shown in Fig- 
ure [T5J In both Newtonian and relativistic cases, the form of the wave is not changed. It 
indicates that the computations are successfully performed in our code. The propagated 
waveforms are very smo oth and no oscillations are found like the ones in the previous study 
JPe Villiers et alil2003h . 
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Fig. 12. — Top panels: Newtonian Alfven wave propagation in z direction. The left figure 
show initial and final profile of v z and the right one does that of B z . Bottom panels: 
Relativistic Alfven wave propagation in z direction. The left figure show initial and final 
profile of v z and the right one does that of B z . 
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C. Conservative Variable and Fundamental Variables 

In our numerical code, variables such as hWv 1 , D, E, B l are evolved. These variables 
are called conservative variables. It is necessary to compute fundamental variables such as 
v' t ,p,e,B l from these conservative variables. In our computations, pressure is not described 
as an analytic function of energy e, therefore algorithms used in other GRMHD simulations 



(IDe Villiers et al.l 120031 ) is not available here. If the value of the enthalpy, h, is found, 
the Lorentz factor, W, is obtained from the values of hWv 1 and then all the fundamental 
variables are obtained. To determine the value of h, we search the root of the equation 
below: 

f(h) = (1 + e/p + p/p +\b\ 2 /p)-h = (CI) 

where \b\ 2 = (B 2 + bfy/W 2 . The fundamental variables in the equation are obtained when h 
is assumed, as stated above. We use simple bisection method to search the root of Eq. IC1I 
and make the error of the equation, Ah/h, below 10~ 4 . 



D. Test Problems 

At last we show some results of the numerical tests on our code. We have done three 
typical problems in both non-relativistic case and relativistic case. Three problems are the 
shock tube problem, reflection shock problem and magnetic shock tube problem. We present 
the results one by one. 



Shock Tubes We consider the shock tube used by Sod (ISodlll978l ) in his comparison of 
finite difference sch eme. First we perfor m weak shock problem and strong shock problem 
in Newtonian case (IHawley et al.l Il984bl ). For the weak one, the initial conditions of this 
problem gas with T = 1.4 with pressure and density Pi = 1.0, pi = 10 5 to the left of x = 0.5 
and P r = 0.1, p r = 0.125 x 10 5 to the right. For the strong one, the initial conditions of this 
problem are gas with T = 1.4 with pressure and density Pi = 0.67, pi = 100 to the left of 
x = 0.5 and P r = 0.67 x 10~ 7 , p r = 1.0 to the right. The numerical values greatly correspond 
to the analytic value. 

Next we perform relativistic shock tube problems. We fix the hyd rodynamical variable s 
except for the pressure of the left side, pi = 10, p r = l,p r = 10~ 6 (lAnninos et al.ll2003l ). 
We set three types of pressure, i.e. pi = 1.33(W = 1.08) for lowly relativistic case, pi = 
6.67(W = 1.28) fir mildly relativistic case and pi = 666. 7(W = 3.28) for highly relativistic 
case. The results are shown in the left panels of Figure [121 For lowly and mildly relativistic 
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case, the numerical value greatly correspond to the analytic value. For highly relativistic 
case, the velocity of the numerical comparison doesn't reach that of the analytic solution. It 
is due to the artificial viscosity which converted kinetic energy to internal energy. 



Wall Reflections A second test presented here is the wall shock problem involving the 
shock heating of cold fluid hitting a wall at the left boundary. When the fluid hits the wall 
a shock forms and travels to the right, separating the pre-shocked state composed of the 
initial data and the post-shocked state with solution in the wall frame 



V s 

f>2 



P2 (T-1)(W 1 -1) 

T + i r 

pi ' 



r- 1 r- 1 



(Wi - 1) 



(Dl) 
(D2) 
(D3) 



where Vs is the velocity of the shock front, and the pre-shocked energy and post-shocked 
velocity were both assumed negligible (e = V2 = 0). 

The initial data are set up to be uniform across the grid with adiabatic index T = 4/3, 
pre-shocked density p\ = 1, and pre-shocked pressure Pi = 10 -6 . And we change the velocity 
of the unshocked region parametrically. The result is shown in the right panels of Figure 
[TBI For all computations, the differences between the numerical solution and the analytic 
one are small, however in case of relativistic one, pressure of the numerical solution is bigger 
than the analytic one. It is also due to the artificial viscosity assumed here. 
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Fig. 13. — Left and right panel show the shock tube and the wall reflection tests, respectively. 
From top to bottom, the Lorentz factor becomes larger (for right panels, top: v = 0.1c 
(W = 1.01), middle: v = 0.5c{W = 1.33), and bottoms = 0.9c(W = 5.26)). Note that the 
pressure and density is normalized by 100. 
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Mag netic Shock Tubes At last we present magnetic shock tube problems (jBrio &: Wu 
19881 ). We show the initial condition and the results of the computations in Table HI And we 
present mildly relativistic case in Figure dU showing that our code can handle the various 
magneto-sonic waves as good as the code by lDe Villiers et al.l (120031 ). 




Fig. 14. — Mildly relativistic magnetic shock tubes 
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Table 4: Initial and intermediate states for shock tube tests. This table list the measured 
values in each state of shock tube: "Left" is the initial left state for the given variables, "FR" 
is the value of the variable at the foot of leftmost fast rarefaction fan, "SC" is the value of 
the peak of the slow compound wave, "CDi" is the value of left of the constant discontinuity, 
"CD r " is the value of right of the constant discontinuity, "FR" is the value of the variable at 
the foot of the second fast rarefaction fan and "Right" is the initial right state for the given 
variables, 



Test 


Variable 


Left 


FR 


SC 


CD] 


CD r 


FR 


Right 


Newtonian 


P 


1.00 


0.66 


0.84 


0.70 


0.25 


0.12 


0.13 


^(xlO- 2 ) 


P(xl0~ 4 ) 


1.00 


0.44 


0.73 


0.50 


0.50 


0.09 


0.10 


= 0.75 


^(xlO- 3 ) 


0.00 


6.67 


4.6 


5.99 


6.02 


-2.79 


0.00 




^(xl0~ 2 ) 


0.00 


-0.25 


-1.10 


-1.58 


-1.58 


-0.20 


0.00 




-zRxlO- 2 ) 

V47T V 


1.00 


0.6 


-0.5 


-0.03 


-0.54 


-0.9 


-1.00 


Mildly 


P 


1.00 


0.65 


0.84 


0.70 


0.24 


0.11 


0.13 


Relativistic 


P(xl0~ 2 ) 


1.00 


0.42 


0.71 


0.49 


0.49 


0.08 


0.10 


^(xlO- 1 ) 


^(xlO- 1 ) 


0.00 


0.67 


0.46 


0.58 


0.58 


-0.26 


0.00 


= 0.75 


^(xlO^ 1 ) 


0.00 


-0.24 


-0.94 


-1.5 


-1.5 


-1.9 


0.00 






1.00 


0.56 


0.3 


-0.52 


-0.52 


-0.88 


-1.00 


Relativistic 


P 


1.00 


0.59 


0.70 


0.65 


0.31 


0.11 


0.13 




P 


1.00 


0.51 


0.60 


0.47 


0.47 


0.08 


0.10 


= 0.75 


v x 


0.00 


0.41 


0.27 


0.28 


0.28 


-0.12 


0.00 




v y 


0.00 


-0.07 


-0.62 


-0.58 


-0.11 


-0.11 


0.00 




By 


1.00 


0.61 


0.19 


-0.24 


-0.24 


-0.81 


-1.00 



